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Abstract 



We calculate the surface plasmon dispersion relations for a periodic chain 
of spherical metallic nanoparticles in an isotropic host, including all multipole 
modes in a generalized tight-binding approach. For sufficiently small particles 
{kd <C 1, where k is the wave vector and d is the interparticle separation), 
the calculation is exact. The lowest bands differ only slightly from previous 
point-dipole calculations provided the particle radius a ;$ d/3, but differ sub- 
stantially at smaller separation. We also calculate the dispersion relations for 
many higher bands, and estimate the group velocity v g and the exponential 
decay length £d for energy propagation for the lowest two bands due to single- 
grain damping. For a/d = 0.33, the result for £d is in qualitative agreement 
with experiments on gold nanoparticle chains, while for larger a/d, such as 
a/d = 0.45, v g and £d are expected to be strongly /c-dependent because of the 
multipole corrections. When a/d ~ 1/2, we predict novel percolation effects 
in the spectrum, and find surprising symmetry in the plasmon band structure. 
Finally, we reformulate the band structure equations for a Drude metal in the 
time domain, and suggest how to include localized driving electric fields in 
the equations of motion. 
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I. INTRODUCTION 



Recently, it has been shown that energy can be transmitted along a one-dimensional 
(ID) chain of equally spaced metallic nanoparticles, via propagating surface plasmon (SP) 
modes [1-8]. These modes are basically linear combinations of single grain SP modes, 
i. e. oscillations of electronic charge within a single grain [9]. The single-grain modes are 
accompanied by an oscillating electric moment (dipole and higher) on the grain. The electric 
field of this moment in turn generates oscillating moments on neighboring spheres. 

The propagating SP modes are simply traveling waves of these oscillating moments. They 
are characterized by well defined dispersion relations u>(k) which relate their frequencies cu 
and wave vectors k [2,3]. If the damping is sufficiently small, the energy transmitted by 
these ways may travel at speeds up to 0.1c. Thus, one can imagine a variety of possible uses 
for these waves. 

The calculation of u(k) typically involves several approximations. The first of these 
is the near field approximation - that is, one assumes that M < 1, where k is the wave 
number and d the interparticle separation. This assumption permits the electric field E to 
be calculated in the quasistatic limit, in which E is expressed as the gradient of a scalar 
potential. 

A second common approximation is that the field produced by a given particle at its 
neighbors is that of a point dipole. However, this second assumption is stronger than the 
quasistatic approximation, and becomes inaccurate when the particles are closely spaced. 
Under these conditions, the quasistatic fields may be modified significantly by multipolar 
interactions. Typically, these multipolar fields have been included by numerical techniques 
such as finite-difference time-domain calculations. However, it may be difficult to obtain 
accurate results by these numerical techniques, because the higher multipole excitations 
may produce fields which vary rapidly in space, whereas numerical studies with insufficient 
discretization may not achieve adequate resolution. Thus, an exact calculation in a sim- 
ple geometry may not only provide physical insights into this system, but also give useful 
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guidelines for the validity of numerical calculations in more complex geometries. 

In the present work, we will show how these multipolar corrections can be calculated 
exactly, using a straightforward analytical approach. The formalism is analogous to the gen- 
eralized tight-binding method in conventional band theory. In this approach, one constructs 
Bloch states from individual atomic orbitals, and then diagonalizes the Hamiltonian matrix 
in the basis of these atomic orbitals. In the surface plasmon analog, the individual atomic 
orbitals are multipolar SP oscillations for each sphere. The matrix elements needed to calcu- 
late the Hamiltonian matrix are easily constructed, especially for a periodic one-dimensional 
(ID) chain of spheres. The diagonalization needed to calculate the bands is readily carried 
out. The entire calculation is made simpler in ID systems, because the Hamiltonian matrix 
decomposes into separate blocks, one for each azimuthal quantum number m. 

The basic formalism necessary to carry out this calculation [10] has thus far been applied 
only rather sparingly, because there have been few experimentally available realizations of 
the ordered geometry required for this approach. It has been applied mainly to calculat- 
ing the effective complex dielectric function e e {uj) of a periodic composite medium, which 
requires the SP band structure only at Bloch vector k = 0. Only recently has it become 
possible to produce well-controlled ordered metallic structures at the nanoscale, and hence 
to generate and detect these SP waves at general wave vectors. In this paper, we describe 
and numerically solve the equation necessary to calculate this band structure in the general 
case of k 7^ in ID systems. 

Besides giving the solutions in the full multipolar case, we also include damping of the SP 
modes due to losses within the individual metallic particles. For a Drude metal, these losses 
can be treated by including a finite relaxation time r in the Drude dielectric function. This 
damping can be included approximately by adding an imaginary part to the SP frequencies. 
We also briefly discuss why the damping due to radiative energy losses is expected to be 
small. These losses arise from the breakdown of the quasistatic approximation, and can, 
in principle, also be approximately included by adding an imaginary part to the surface 
plasmon frequency. 
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We will also present the multipolar SP equations in the time domain for a Drude metal, 
where they take the form of a set of coupled second order ordinary differential equations. 
In this form, it is straightforward to include single-particle damping (and also, in principle, 
radiative damping). Moreover, one can also incorporate driving terms, arising, e. g., from 
external electric fields. These equations may thus be useful in modeling specific types of 
experimental probes which produce localized time-dependent electric fields. 

The remainder of this paper is organized as follows. In Section II, we describe the 
formalism needed to calculate the SP band structure, and specialize to the calculation for 
a ID chain. Section III presents numerical results in this ID chain system. In Section IV, 
we discuss our results, and suggest some possible extensions to other geometries. Finally, 
the Appendix presents an alternative formulation of the equations of motion in the time 
domain, and describes how localized time-dependent driving electric fields can be included 
in these equations. 



II. FORMALISM 

We consider a ID chain of spherical particles of radius a, separated by a distance d 
(d > 2a). The particles are assumed to be arranged along the z axis with centers at z — 0, 
±d, ±2d, ■ ■ ■. We assume that the particles and host have dielectric functions e m (oo) and 
€h(ou). To be definite, we may consider the particles as metallic and the host as insulating, 
but the discussion below applies to any choice of e m and e h . All the formalism is given in 
terms of a frequency variable s defined by 

s= (1) 

As will be seen below, and as is discussed indirectly in Ref. [10], all the allowed propagating 
SP frequencies correspond to s in the range < s < 1, or equivalently, to — oo < e m /e h < 0, 
assuming that e m and are both real. 

We will calculate the SP band structure in an "atomic" basis n, £, m. Here, £ and m label 
the "angular momentum" of the eigenfunction, and n labels the grain. Thus the allowed 
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values of these indices are n — 0, ±1, ±2, • • •, £ = 1, 2, 3, • • •, and m = —I, — £ + 1, • • • , I. 
In order to calculate the SP band structure in this basis, we first need the matrix element 
Qnfrn;n't'm' , where n ^ n' . This matrix element is given by (see, e. g., Ref. [10]) 



Qntm;n'l'rn' ~ ( _ 1) 



t'+m' 



e+e'+i i v 1/2 



\n-ri\d) V( 2 ^ + 1)(2£' + 1), 



^(^ +TO )i(iVmV(r-^- TO oi] i / a - exp(#6K ~ m))p ^ m(c ° s ^ (2) 

Here we have introduced b = {p! — n)dz, which is the vector separation between the grains 
centered at n'z and nz, and the polar and azimuthal angles 9b and (f>b for this displacement 
vector. Since b = bz, 9b is either or n, depending on whether b is positive or negative. 
If b > 0, P™+7 m = 5 m ', m , whereas if b < 0, P™^ m = (-l) e ' +e 5 m ', m - Incorporating these 
simplifications, we find that with 

/ „ \ W+ 1 / fifii \ V2 



Qnim-.n'i'm' ~ (~ 1 



(£ + f )! ( n' -n V +e ' 

X p + mjUj' + m)\(£- mjUji - m)!] V2 [ \ n > - n \) 6m > m ^ 

The matrix elements are diagonal in m because the one-dimensional chain is unchanged on 
rotation by any angle about the z axis. 
Next, we define the matrix element 

Qtm;i'm'{k) = Qo,£m;nfJm exp(inkd)5 m ,m> , (4) 

where the sum runs over all positive and negative integers except n — 0. We have used the 
fact that, for this periodic ID system, Q n em;n'e'm' is a function only of n — n' and vanishes 
for m 7^ m'. After a little algebra, using Eqs. (3) and (4), we obtain 

n (h\ f*\ w ' +1 ^ ™<nkd) T ^ 

Wim;i'm'\K) — I - I 2^ I+F+l ,mO m ,m> [O) 

\aj n=1 n 



for I + £' even, and 



n (b\ fa\ W ' +1 ^ Mnkd) T ^ 



for £ + £' odd, where 

K 2( i ( te \ V2 + m 

e ' i '* m y ' \(2£+l)(2£> + 1)) [(£ + m)\(£> + m)\(£-m)\(£> -m)!.} 1 / 2 ' K ) 

Finally, in terms of the matrix elements Qemi'm] the SP structure is given by 

det\s-H\=0, (8) 
where the "Hamiltonian" H has matrix elements 

He m ;e'm(k) = Se5e,e> + Qtm;i>m{k), (9) 

and the "atomic" eigenvalues se are given by 

s ' = 2fTi- (10) 

Note that s 1 = 1/3, while s e — > 1/2 as £ — > oo. 

The full SP band structure at wave vector is obtained by diagonalizing the matrix (9). 
As in conventional band theory, there are many band energies for a given k, and one need 
consider only the bands in the first Brillouin zone, i. e., in this case, for — n/d < k < n/d, 
since the states at other values of k are equivalent to those in the first zone. For this one- 
dimensional system, the Hamiltonian breaks into separate blocks, one for each value of m; 
this conveniently reduces the size of the matrix which needs to be diagonalized. Finally, as 
in the LCAO (linear combination of atomic orbital) method of conventional band theory, 
the band structure that results from this analysis is composed of bands which originate from 
various "atomic" orbitals. In the present case, the the "atomic" states are multipolar SP 
modes associated with the individual spheres. These are degenerate at different m (since 
the individual particles are spheres), and have eigenvalues S£ = £/(2£+ 1). 

The band structure that results from diagonalizing the matrix (9) is expressed in terms 
of the variable s. Thus, the bands have the form s a (k), where a labels the individual 
bands. These may be converted into frequencies using the relation (1). This dispersion 
relation will take on various forms, according to how e m (and e h ) depend on uj. Here we 
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assume that the system consists of Drude metal particles in vacuum, so that = 1 and 
e m {uj) = 1 — u)p/[u{uj + i/r)]. If u p t — > oo, then the appropriate conversion is given by 

u a (k) = u p \]s a {k). (11) 

For £ = 1 there are three degenerate modes at frequency uj p / \/3, while for t — > oo, the modes 
approach the limiting value of tu p / \/2. 

If the relaxation time is finite, then the relation 

, = = i^L±iA) (12) 

l-e m /e h 

can be inverted to express u as a function of s. For a real s, the corresponding uj has both a 
real and imaginary part. Thus, we can map an eigenvalue s a {k) (where a is a band index) 
to a complex eigenvalue 

u a (k) = ^ p s a (k)-l/(Ar 2 )-i/(2r). (13) 

The imaginary part describes the damping of this mode due to the finite lifetime of the 
surface plasmons in the individual spheres. If u p r ^> 1, this damping is small, and the 
shift due to the damping (as embodied in the factor 1/r 2 within the square root) is even 
smaller. Note that we have not included radiative damping in this expression. In contrast 
to single-particle damping, the radiative damping depends on the particle size, being greater 
for larger particles. For 10-nm radius gold spheres, it has been estimated that only 1.5% of 
the total damping rate is due to radiative damping [11]. Also, according to Refs. [3] and [4], 
the radiative damping should be small for particles of such a size because of strong near-field 
interactions. 

In the Appendix, we present an alternative formulation of the equations of motion in the 
time domain to obtain the SP band structure, assuming a Drude dielectric function. In this 
formulation, which has the advantage that it can deal with localized time-dependent driving 
electric fields, these radiative corrections could easily be included, as has been discussed, for 
example, in Ref. [3]. 
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In order to compare with experiment, we will consider two more quantities which can be 
obtained from the dispersion relation u>(k): First, the /c-dependent group velocity v g (k) is 
given by the relation 

vM = % (14) 

and can be easily computed numerically, given u(k). Secondly, we can also use Eqs. (13) 
and (14) to estimate the energy loss from a plasmon propagating along a chain, which is 
important in applications. For this purpose, we define energy decay lengths, £d(&) for the 
lowest longitudinal and transverse modes, as the distance over which the energy density in 
the wave amplitude decreases by a factor exp(— 1). If the complex band frequency is denoted 
u>i(k) + iu 2 (k), then is defined by 

_ v g (k) _ 1 du(k) 
d 2uo 2 d 2uo 2 d(kd)' 1 ' 

Note that in the case of the Drude approximation, the imaginary part of the complex band 
frequency does not depend on k, and thus £d is just proportional to v g (k). 



III. NUMERICAL RESULTS 

We have diagonalized the matrix (9) to obtain the surface plasmon band structure for 
various values of the parameter a/d. We include all bands up to i = 80, which is sufficient to 
insure convergence of s a (k) to within 1%. The results are shown in Fig. 1. For comparison, 
we also show the results of truncating the Hamiltonian matrix at £ = 1. In this latter case, 
the Hamiltonian matrix is a diagonal 3x3 matrix with elements 

/,\ 1 f a \ 3 cosinkd) T ^ 
H lm ., lm (k) = 3 + (j) E ^3 X i,m, ( 16 ) 

where i^i,i,o = —4/3 and i^i,i,±i = 2/3. The corresponding plasmon bands, expressed as 
s a (k), are shown as open square (m = 0) and circle (m = ±1) in Fig. 1. If we use the 
Drude expressions w^(k) = UpS a (k), then these correspond exactly to the dipolar SP band 
structures obtained in Ref. [3]. This behavior is as expected, since when we retain only the 
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£ = 1 terms in the Hamiltonian, we are neglecting all the quasistatic contributions except 
the dipole fields. 

As is evident, the plasmon bands take on quite a different form when the higher values of 
£ are included in the Hamiltonian matrix. For small values of a/ d (i. e. a/d < 0.3), the lowest 
three bands are very similar to the bands obtained when only the £ = 1 matrix elements are 
included. This behavior is not surprising, since at these values of a/d the matrix elements 
connecting the £ = 1 states to those of higher £ are quite small, for any value of k. This 
smallness originates in the factors of (a/d) i+i ' +1 which appear in all the expressions for the 
matrix elements. The smallest factor connecting £ = 1 to higher £ is (a/d) 4 for k ^ 0, and 
(a/d) 5 for k — 0. Thus, for relatively small values of a/d, these connecting matrix elements 
are substantially smaller than the diagonal ones. 

When a/d > 0.35, the interband matrix elements start to become substantial. When 
this happens, the shapes of the lowest bands start to depart significantly from the purely 
dipolar form seen for smaller a/d. As is evident, by the time a/d — > 1/2, the band structures 
of the lowest bands are so altered that they no longer bear any obvious relation to these 
dipolar band shapes. Precisely at a/d = 1/2, the lowest state at k — reaches the limiting 
value s — 0. This behavior is a percolation effect: when a/d = 1/2, the two neighboring 
spheres just touch, and the entire line of spheres becomes one connected chain. Thus, one 
might expect that the lowest eigenvalues of this chain would resemble that of a very long 
cylinder, which indeed has as its lowest eigenvalue s — 0. 

The band structure also acquires a striking symmetry near a/d = 1/2. First, there 
appears to be a nearly perfect reflection symmetry about the line s = 1/2. In addition, 
there is another reflection symmetry about the line k = n/(2d), i. e. at the middle value of 
k in the first Brillouin zone. As particular examples of these symmetries, there appear to 
be eigenvalues of s = and s = 1 at k = n /d, just as there are at k — 0. We do not fully 
understand the reasons for these symmetries. The s ~ eigenvalue at k — ir/d apparently 
corresponds to a longitudinal mode (dipole moment of the spheres parallel to the z axis) in 
which each sphere oscillates 180° out of phase with its neighbors. The multitude of modes 
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near s = 1/2 presumably originate in the high-£ "atomic" modes, which have eigenvalues 
approaching s = 1/2. 

In Fig. 2, we show the eigenvalues of the two lowest bands at k = plotted as a function of 
a/d. Here, the lowest band corresponds to longitudinal mode (m = 0) and the second lowest 
band to degenerate transverse modes (m — ±1). We performed two different calculations: In 
the first calculation, shown as open circles and squares, we assumed the dipole approximation 
and included only the £ = 1 part of the Hamiltonian matrix. This calculation corresponds 
to the tight-binding approximation used in Ref. [3]. In the second, we included all bands up 
to £ = 80, which is sufficient to ensure the convergence of these two bands, as in Fig. 1, and 
this inclusion of the higher multipoles starts to make a substantial difference for a/d> 0.35. 

The inset to Fig. 2 shows the splitting As between the longitudinal and transverse 
modes at k — 0, plotted as a function of a/d. In the dipole approximation (dashed line), 
this splitting increases monotonically as a/d increases. However, as shown by the solid line 
in the inset, when the higher multipoles are included, the splitting reaches a maximum near 
a/d = 0.46, then decreases again. 

In order to compare with experiment, one needs to re-express the band structure as 
dispersion relations for uo(k) rather than s(k) ) using Eq. (1). With the resulting uj(k), we 
can also obtain v g (k) from Eq. (14) and £d from Eq. (15). We show the resulting dispersion 
relations u(k)/u p in Fig. 3 (a), and the resulting £d(/c) and v g (k) in Fig. 3 (b) for the lowest 
longitudinal and transverse bands as a function of kd. We denote these results as open 
square and circle for a/d = 0.33, the value used in experiments, and also as solid and dashed 
lines for a/d = 0.45, which is near the maximum of the splitting As. In order to calculate 
u(k)/ujp, we choose u p = 6.79 x 10 15 rad/s and r = 4 fs, as used in Ref. [7]. This choice 
allows us to compare the present result for a/d = 0.33 with those given in Refs. [5] and [7]. 

First, we compare our results for a/d = 0.33 with experiment. For a/d = 0.33, the result 
of the full calculation for u(k) is not significantly different from the dipole approximation, 
since multipole effects produce only a minor alteration to the lowest bands in this case. 
However, the multipole effects can be seen much more clearly in the /c-dependent group 
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velocity v g (k) for these bands, and this quantity can be easily computed numerically using 
the u){k) shown in Fig. 3 (a). In contrast to the result from the dipole approximation of Ref. 
[3], the maximum in v g for a/d = 0.33 does not occur at k — n/(2d), but instead around 
k = it /(Ad), when the multipolar corrections are included. However, if we assume d = 75 
nm, Up = 6.79 x 10 15 rad/s, and r = 4 fs as in Ref. [7], the magnitude of the maximum v g for 
the longitudinal (m = 0) mode is approximately 1.9 x 10 7 m/s, which is close to the result 
of Ref. [7], while the magnitude of v g for the transverse (m = ±1) modes is slightly larger 
than the value (1.1 x 10 7 m/s) estimated in Ref. [7]. Also, the values of £d in the lowest 
longitudinal and transverse modes for a/d = 0.33 are comparable to the experimental values 
for gold, as given in Ref. [6]. 

For a/d = 0.45, which is near the maximum of the splitting As, the multipole corrections 
to the band structure are much greater. As a/d approaches the maximum splitting between 
longitudinal (m = 0) and transverse (m = ±1) modes, the variation of v g with k becomes 
non-monotonic. In contrast to the dipole approximation, which gives a maximum group 
velocity at k — n/(2d), our exact calculation actually gives a local minimum in v g for this 
value of k (for both polarizations). As can be seen from Fig. 3(b), v g (k) has, in fact, two 
maxima as a function of k for this separation, for both longitudinal and transverse modes. 
The maximum estimated exponential decay length shown in Fig. 3(b), for the optimum k, 
corresponds to an m = wave, and is about three times larger than that for a/d = 0.33. 
But this decay length is calculated for a wave with k vector corresponding to the maximum 
group velocity. The actual v g is strongly /c-dependent, especially for the larger a/d. Thus a 
typical wave, which would likely propagate as a packet of many different wave vectors, would 
likely have a quite different decay length, and also would probably not decay exponentially. 
It is possible that this /c-dependence is related to the non-exponential spatial decay of SP's 
found in the numerical simulations of Ref. [1]. 

We have not commented thus far about the role of the higher SP bands. For values of 
a/d greater than about 0.33, most of these bands are nearly dispersionless, with eigenvalues 
s a (k) ~ 1/2. The SP modes corresponding to these bands will thus propagate with very 
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small group velocity v 9t0l = du a (k)/dk, and are likely to contribute very little to energy 
transport along the chain. 



IV. DISCUSSION 



In this work, we have calculated the surface plasmon dispersion relations for a one- 
dimensional chain of spherical particles in a uniform host. In contrast to previous work, we 
have included all the multipolar terms in our calculation, within the quasistatic approxi- 
mation. We find that the dipole approximation is reasonably accurate for a/d < 1/3, but 
becomes increasingly inaccurate for larger a/d. When a/d — > 1/2, the lowest band shapes 
are entirely different from the point-dipole predictions. Thus, an accurate comparison of 
theory to experiment should take into account these corrections, when a/d exceeds about 
1/3. 

In our results near a/d = 1/2 we see conspicuous percolation effects. Specifically, the 
k = mode approaches zero frequency for a chain of Drude spheres in an insulating host. 
Furthermore, when a/d — > 1/2, the entire band structure shows remarkable reflection sym- 
metry, both around k = 7r/(2d) and around the frequency midpoint at s = 1/2. We do not 
presently understand the reasons for this symmetry. 

Besides producing shape distortions in the lowest bands, the present calculations also 
lead to an infinite number of higher propagating SP bands. We believe, however, that these 
will contribute little to energy propagation, because they are characterized by much lower 
group velocities than the lowest two bands. 

Our calculations in Fig. 3 have, of course, been carried out in the Drude approximation. 
As mentioned earlier, we used values of 1/r and u p as best fits to experiments on bulk gold, 
as described in Ref. [7]. In actuality, the complex dielectric functions of silver, and especially 
gold, have substantial interband contributions, and cannot be described by a Drude dielectric 
function in the visible. An accurate translation of the SP band structure from the variable s 
to the variable e m (cj)/e/ l should use this more accurate dielectric function, e. g. by using a fit 
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of the experimental e m (u) to a sum of free-electron and Lorentz oscillator parts. This more 
complicated procedure might somewhat change the plasmon band structures, especially for 
gold. Another possible complication is that, in typical experiments, the nanoparticle chains 
are laid down on a substrate, whose dielectric constant differs from that of vacuum. Thus, 
the chain is not really embedded in a homogeneous dielectric. Some workers have taken 
this complication into account by treating the host as homogeneous but with a dielectric 
function which is an average over the air and substrate dielectric functions [5]. Once again, 
this correction, if included, would also modify the calculated SP band structure. 

The present work could be readily be generalized to higher dimensions, e. g. to an 
ordered layer of spheres deposited on a substrate. This extension should be straightforward, 
since the matrix elements required are the same as used here. The same approach could 
also be applied to particles of different shapes, e. g. ellipsoids or short cylinders, although 
the calculation of the single particle eigenstates and the overlap integrals might be more 
difficult. We plan to carry out some of these extensions in the near future. 
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APPENDIX: FORMULATION IN THE TIME DOMAIN 

If the small spherical particles are described by a Drude dielectric function, the SP band 
structure can be also obtained, in perhaps a more intuitive way, by writing down a set of 
coupled equations of motion in time for the multipole moments. We first introduce the 
scalar quantities q n e m , defined as the (£m) th multipole moment of the n th particle. Then, in 
the absence of damping, the coupled equations of motion can be written in the form 

Qnim — —UfmQnem + Q Vm;n' '£> 'm> 'Qn' '£> 'm' (17) 

e'm'n' 
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where the prime over the sum indicates a sum only over the grains n' ^ n. For spherical 
grains, the single-grain resonant frequencies are given by 

and the coupling constants Qntmyi'm 1 are given by Eq. (3). Eqs. (17) are readily solved for 
the eigenfrequencies by substituting assumed solutions of the form 

q n tm(t) =Re[qe m exp(inkd - iuj a t)\ (19) 

into Eqs. (17). Here qt m is the amplitude of the {tm) th multipole in a propagating mode 
of wave vector k. With this substitition, Eqs. (17) reduces to a set of coupled homoge- 
neous linear algebraic equations. A solution is obtained if the determinant of the matrix of 
coefficients vanishes. This condition is equivalent to that of Eqs. (8), (9), and (11). 

Eq. (17) has a straightforward physical interpretation. The right-hand side of Eq. (17) 
describes two contributions to the force acting on the multipole moment q n £ m . The first 
term is the restoring force due to charge motion within a single particle. The second term 
on the right comes from the electric fields of all the multipole fields from the other particles, 
evaluated at the position of the n th particle. Damping is easily included in Eq. (17) by 
adding to the right-hand side the term —q n em/T. To obtain solutions in the presence of 
damping, we assume the form (19) but with a complex frequency cu a (k) = uj\ + iu 2 . The 
resulting uj a {k) is given by Eq. (13). 

An appealing feature of Eqs. (17) is that one can easily add a driving term. For example, 
if a uniform electric field E n (t) is applied to the n th grain, the interaction energy between 
that field and the n th grain would be H' = — p n • E n (t), where p n is the dipole moment of 
the n th grain. To calculate the force on the n th grain, we write p n = gx n , where q is the total 
electronic charge in the n th grain, and x n is its displacement from its equilibrium position. 
The interaction energy between this charge and the applied field is thus — gx„ • E n (t). The 
corresponding force on the charge is just gE = Mx n , where M is the total mass of the 
electronic charge in the grain. Thus p n = gx n = jj^E = (47ra 3 n e e 2 /3m e )E = (a 3 ^p/3)E, 
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where m e = 3M/(47ra 3 n e ) is the electron mass, e = 3g/(47ra 3 n e ) is the magnitude of the 
electron charge, n e is the electron density, and u 2 = 4nn e e 2 /m e is the squared plasma 
frequency. 

Finally, to incorporate the damping and driving terms into the right-hand side of Eq. 
(17), we express the applied electric field in terms of the spherical components m = and 
m = ±1. Thus, we write E n ^ Q (t) = E njZ (t), E nAj±1 (t) = E n>x (t) ± E n>y (t). We then obtain 
the following equations of motion: 

1 Up" ' CC^ i 

(jnim = —^tmQntm Qnim H — E n ^ m {t)8^i + UJ p Q nlm;n> I 1 m' In' i> m' ■ (20) 

T 6 n'e'm' 

Eqs. (20) are generalizations of the equations written down in Ref. [3] which include all 
multipole moments, within the quasistatic approximation, and single-grain damping within 
the Drude approximation. They also include the effects of a uniform but time- dependent 
electric field applied to the n th nanoparticle. 

Finally, we note that we have not included radiative corrections to the calculated SP 
bands. However, the present work could also be extended beyond the quasistatic approxi- 
mation to include radiative corrections, in a simple manner, by adding an additional imag- 
inary part to the eigenvalues. These corrections could easily be included within the dipole 
approximation, as has been discussed, for example, in Ref. [3]. 
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FIGURES 




k k 

FIG. 1. Dispersion relations s(k) for the surface plasmon bands propagating along a chain of 
spherical nanoparticles of dielectric function e m in a host of dielectric function eh, plotted versus 
wave vector k. (a) a/d = 0.25, (b) a/d = 0.33, (c) a/d = 0.4, (d) a/d = 0.45, (e) a/d = 0.49, 
(f) a/d = 0.5. The solid and dashed curves correspond to m = and m = ±1 respectively for 
the full band structure, incorporating all bands up to £ = 80, as obtained diagonalizing the full 
Hamiltonian matrix [Eq. (9)]. The open squares (m = 0) and circles (m = ±1) denote calculations 
for the t = 1 modes in the dipole approximation. Note that the m = ±1 modes are degenerate. 
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FIG. 2. Eigenvalues s(k) for the lowest two bands of the band structure shown in Fig. 1 at 
k = 0, evaluated as a function of a/d. Solid and dashed curves correspond to m = and m = ±1 
respectively for full multipole calculations. Open squares (m = 0) and circles (m = ±1) are 
point-dipole calculations. Inset: Splitting As between lowest m = ±1 and m = bands at k = 
as calculated in the dipole approximation (open squares) and using full Hamiltonian matrix (solid 
line). 
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FIG. 3. (a) Dispersion relations u>(k) for the lowest two bands in a chain of metallic nanoparti- 
cles at a/d = 0.33 or 0.45. The solid (m = 0) and dashed (m = ±1) lines correspond to a/d = 0.45; 
the open squares (m = 0) and circles (m = ±1), to a/d = 0.33. The curves are computed using 
the full Hamiltonian up to £ = 80, using a Drude dielectric function for the metal, (b) Energy 
decay lengths £d, in units of the lattice constant d, and corresponding group velocities v g in units 
of ujpd, plotted versus kd for the lowest two bands, assuming a/d = 0.33 or 0.45. The labeling of 
the curves follows the notation of Figs. 3 (a). 
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